function phi_n = cfdInterpolateFromElementsToNodes(phi)
%--------------------------------------------------------------------------
%
%  Written by the CFD Group @ AUB, Fall 2018
%  Contact us at: cfd@aub.edu.lb
%==========================================================================
% Routine Description:
%   This function interpolates the field from elements to nodes
%--------------------------------------------------------------------------

% Get field type
theNumberOfComponents = size(phi,2);

% Get info
theNumberOfInteriorFaces = cfdGetNumberOfInteriorFaces;
theNumberOfElements = cfdGetNumberOfElements;
theNumberOfNodes = cfdGetNumberOfNodes;

theNodeCentroids = cfdGetNodeCentroids;
theNodeFaceIndices = cfdGetNodeFaceIndices;% 围着nodes的面的全局索引
theNodeElementIndices = cfdGetNodeElementIndices;% 围着nodes的elements的全局索引

theElementCentroids = cfdGetCentroidsForElements;
theFaceCentroids = cfdGetFaceCentroidsSubArrayForFaces;

% Initialize node array
phi_n = zeros(theNumberOfNodes, theNumberOfComponents);

% Loop over nodes
for iNode=1:theNumberOfNodes
    
    % Current node properties
    N = theNodeCentroids(iNode,:);
    
    nodeFaces = theNodeFaceIndices{iNode};% 围着这个iNode的面的全局索引
    nodeElements = theNodeElementIndices{iNode};% 围着这个iNode的elements的全局索引

    % Initialize
    localPhiNode = zeros(1,theNumberOfComponents);
    localInverseDistanceSum = 0;

    if isempty(nodeFaces(nodeFaces>theNumberOfInteriorFaces))% 围着这个Node的面都不是边界面
        localElementIndices = nodeElements;

        for iElement = localElementIndices % 遍历周围的elements
            C = theElementCentroids(iElement,:);% element的质心坐标

            d = cfdMag(N-C);% element质心和vertex(node)之间的距离
            localPhi = phi(iElement,:);

            localPhiNode = localPhiNode + localPhi/d;% 反距离加权
            localInverseDistanceSum = localInverseDistanceSum + 1/d;       
        end
    else % 围着这个Node的面有部分是边界面
        localBFacesIndices = nodeFaces(nodeFaces>theNumberOfInteriorFaces);% 获取边界面的全局索引
        for iBFace = localBFacesIndices % 遍历所有边界面
            C = theFaceCentroids(iBFace,:);% 边界面的质心
            iBElement = theNumberOfElements+(iBFace-theNumberOfInteriorFaces);% 边界面对应的边界element的全局索引

            d = cfdMag(N-C);% 边界面质心和vertex(node)之间的距离
            localPhi = phi(iBElement,:);

            localPhiNode = localPhiNode + localPhi/d;
            localInverseDistanceSum = localInverseDistanceSum + 1/d;       
        end
        % FIX me: 是不是有遗漏？围着这个Node的面有部分是边界面，同时也有可能有部分非边界面，不用处理了吗？
    end
    phi_n(iNode,:) = localPhiNode/localInverseDistanceSum;
end
